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ABSTRACT 

This paper presents a mathematically rigorous, subspace projection-based reduced order modeling (ROM) 
methodology and an integrated framework to automatically generate reduced order models for spacecraft 
thermal analysis. Two key steps in the reduced order modeling procedure are described: (1) the 
acquisition of a full-scale spacecraft model in the ordinary differential equation (ODE) and differential 
algebraic equation (DAE) form to resolve its dynamic thermal behavior; and (2) the ROM to markedly 
reduce the dimension of the full-scale model. Specifically, proper orthogonal decomposition (POD) in 
conjunction with discrete empirical interpolation method (DEIM) and trajectory piece-wise linear 
(TPWL) methods are developed to address the strong nonlinear thermal effects due to coupled conductive 
and radiative heat transfer in the spacecraft environment. Case studies using NASA-relevant satellite 
models are undertaken to verify the capability and to assess the computational performance of the ROM 
technique in terms of speed-up and error relative to the full-scale model. ROM exhibits excellent 
agreement in spatiotemporal thermal profiles (<0.5% relative error in pertinent time scales) along with 
salient computational acceleration (up to two orders of magnitude speed-up) over the full-scale analysis. 
These findings establish the feasibility of ROM to perform rational and computationally affordable 
thermal analysis, develop reliable thermal control strategies for spacecraft, and greatly reduce the 
development cycle times and costs. 
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1. INTRODUCTION 


Computational models and analysis are increasingly used to support the engineering of spacecrafts at 
every stage of the development [1-4], including design, test, ground-operation simulation, and controller 
development. Every single element in the spacecraft needs to be thermally interrogated to ensure that its 
temperature, thermal distortion and stability remain in the pre-defined operating envelope and in-orbit 
performance requirements are met during all mission phases. Trends in recent years have been towards 
larger thermal models, and have therefore, placed additional computational demands on thermal 
engineers. Attempts to verify designs by modeling and analysis rather than testing further add to this 
burden. 

Currently, there are three categories of models that are prevalently used in spacecraft thermal design: 
general thermal analytical models, approximate thermal analytical models [3, 5, 6], and heuristic models 
[7-9]. In general thermal analytical models, the discretized form of the governing PDEs for heat transfer, 
provide detailed and accurate temperature distributions in the spacecraft that can be used to resolve the 
thermal interaction with other physics for coupled analysis. Concurrently executed user logic and/or other 
customized inputs are typically used as the boundary conditions and heat sources. However, the 
computational process of general models normally is computationally intensive and time-consuming, 
which could take anywhere from days to weeks to yield reasonable answers to design questions. In 
contrast, approximate models based on the lumped-parameter assumption allows fast simulation speed 
and the ability to interrogate critical physics by parametric, tradeoffs, and interface studies at an early 
stage of the thermal design when the spacecraft concept, structure, and test arrangement are yet fully 
defined. Therefore, it provides valuable information about the feasibility of the spacecraft from the 
thermal perspective [6]. However, approximate models are generally poor in accuracy and lack 
quantitative predictive capability. More importantly, existing approximate thermal models are not capable 
of providing detailed temperature distribution in spatial domains; therefore they are inapplicable to 
coupled, multi-disciplinary spacecraft analysis. 

To overcome aforementioned limitations, data-driven heuristic models have been developed for rapid 
parametric analysis and design evaluation. Miller [7, 10, 11] used the response surface methodology to 
study the thermal performance of the Orion vehicle with varying attitudes. The model showed good 
agreement with the general thermal analytical models and quickly screened a large number of verification 
tasks. Prince [9] derived response surface equations from the high-fidelity thermal model of the Mars 
Reconnaissance Orbiter and integrated those equation into the autonomous aerobraking simulation 
software. A goodness of fit analysis was performed confirming the response surface equations adequately 
represented the high-fidelity thermal model. However, issues associated with data driven heuristic models 
include: (1) they are the regression-based method, leading to poor accuracy relative to the mathematically 
rigorous, general thermal models, (2) the regression mapping is established between the input parameters 
and point-wise quantities of interest (e.g., max/min at critical locations), and hence, cannot be used to 
predict the entire spatial temperature distribution within the spacecraft; and (3) in general, they are ill- 
suited for predicting the dynamic thermal behavior of the spacecraft. Therefore, there is a strong need for 
an efficient spacecraft thermal modeling and analysis methodology that enables salient analysis accuracy 
similar to high-fidelity general thermal models but costs a fraction of their simulation time (equivalent to 
the approximate or heuristic models). 

In this context, this paper presents a generalized, mathematically formal, subspace projection-based 
reduced order modeling (ROM) methodology for spacecraft thermal analysis. The computationally 
efficient ROM is automatically derived by projecting a full-scale general thermal model onto a low- 
dimensional, characteristic subspace. Two rigorous algorithms used to construct the subspace, including 
proper orthogonal decomposition (POD) [12-15] and trajectory piecewise linear reduced order modeling 
(TPWL-ROM) [16-19] are presented, and their formulation specific to spacecraft thermal model is 
described. As a result, a low-dimensional ROM whose size is 1-2 orders of magnitude less than the 
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original full-scale model can be obtained. The ROM is then coupled with ordinary differential equation 
(ODE) and differential algebraic solver (DAE) solvers to resolve the spatiotemporal temperature profile 
of the whole spacecraft. The ROM is verified against the full-scale, general thermal models of the 
spacecraft, including the Laser Interferometer Space Antenna (LISA) satellite and the Global 
Precipitation Measurement (GPM) satellite model available at the NASA Goddard Space Flight Center. 

The facets that clearly distinguish the present work from the existing efforts include: (1) the first attempt 
to develop ROMs for spacecraft thermal analysis based on the mathematically rigorous, subspace 
projection method; (2) our ROMs bridge the gap among the existing techniques, and saliently combine 
several benefits of them in terms of computational accuracy and speed. Different from the approximate 
and heuristic model, our ROM is capable of capturing the entire temperature profile within the spacecraft 
in a transient manner, and allows for comprehensive spatiotemporal inspection of the thermal transport at 
the system level. Relative to the full-scale general thermal model, typically 1-2 orders of magnitude 
dimension reduction and computational speedup can be achieved without appreciably compromising 
simulation accuracy. The judicious combination of these merits make ROM-based simulation 
methodology well-suited for initial concept evaluation and screening, orbital analysis, operational 
optimization, as well as development of reliable thermal management strategies; and (3) in addition to the 
ROM algorithms, the present paper also formulates a systematic approach for acquiring and converting 
full-scale models in complex scenarios (e.g., temperature dependent parameters) to the form amenable to 
ROM implementation, laying a firm foundation for subsequent development. 

The paper is organized as follows: the ROM approach and procedure is first introduced in Section 2, in 
which the acquisition of the full scale model and the POD and TPWL algorithms are elucidated. In 
Section 3, the ROMs are verified and demonstrated using two relevant satellite model examples. The 
paper is finally summarized in Section 4. 

2. ROM METHODOLOGY 

In this section, we will first present the systematic organization of the ROM methodology for spacecraft 
thermal analysis. Then the key elements of the full-scale model acquisition and dimension reduction for 
generating reduced thermal models will be described. As shown in Figure 1, the ROM-based thermal 
analysis includes three key steps: 

1. Acquisition of full-scale, general thermal models : The first step is to import the definitive information of 
the spacecraft model from relevant thermal analysis software into the ROM platform via an external 
interface. The information includes conductive and radiative links, node capacitance, environment 
sources, user sources, etc. obtained from appropriate differencing schemes (e.g., finite difference or 
finite element). The model information is then assembled to form a set of ODEs/DAEs in compliance 
with the energy conservation on each node and the connectivity among nodes. The ODE/DAE set hence 
represents the full-scale, general thermal models that can be computed to predict the temperature 
distribution in a transient manner. More details are provided in Section 2.1. 

2. Generation of reduced order model : The assembled full-scale thermal model is then reduced using the 
mathematically rigorous ROM algorithm. Specifically, a low-dimensional subspace is constructed, 
onto which the full-scale model is projected leading to a low-dimension thermal ROM. The ROM can 
be integrated using ODE/DAE solvers to resolve the transient temperature profile in the reduced 
domain. Due to its significantly reduced dimension relative to the full model, the ROM can be 
computed quickly and efficiently. The results of the generated ROM can then be compared against 
the full-scale thermal model for verification. More details are provided in Section 2.2. 

3. Use of ROMs for spacecraft thermal analysis : It should be noted that in many instances, ROM 
generation is a one-time cost, and the validated ROM can be reused in subsequent iterative simulation 
sweeping the operating parameter space to identify the best parameter combination for design 
optimization and thermal control strategy. In contrast to the heuristics-based approach that establishes 
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the point-wise input-output mapping, the projection-based ROM bears excellent accuracy in the 
parameter space to predict the entire dynamic temperature distribution because of its mathematically 
rigorous nature. More details are provided in Section 3. 
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Figure 1. Schematic of the ROM methodology for spacecraft thermal analysis. Model definitive 
information is obtained from the full-scale model via an external interface. The imported model is 
assembled and reduced using mathematically rigorous algorithms to obtain the reduced order model. 
Both full-scale and ROM models are solved and the ROM performance characterized in terms of 
computational speedup and accuracy. Note that due to the use of projection methods the comparison 
between ROM and full-scale model is made over the entire spatiotemporal domain. 


2.1 Acquisition of Full-scale Thermal Models 

The procedure starts with acquisition of the full-scale thermal model. First, a custom interface is used to 
import the definitive information (such as thermal links, heat sources, thermal capacitance, etc.) from the 
geometrical mathematical model provided by the commercially available full-scale thermal analysis tool. 
Next model information is assembled into a dynamic ODE/DAE system governing the heat transfer at each 
individual node (or element depending on the discretization schemes used) in the spacecraft along with 
appropriately selected system inputs 
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( 1 ) 

where 7^, C/, and Qt are, respectively, the temperature, thermal capacitance, and heat source of the node 
(/■=1, 2...n) and n is the total number of the node. Note that Qi includes heat sources from internal heat 
generation (e.g., heaters and electronics heat dissipation) as well as environmental fluxes (e.g., direct solar 
radiation, albedo and earth radiation) [6]; Ky, Ry, are, respectively, the conductive and radiative links 
between nodes. Ni and Mi are, respectively, the number of the neighboring node interacting with the 
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node via conduction and radiation. Ksi and Rsi are the conductive and radiative links between the node 
and an isothermal object (e.g., space), respectively. Ts is the temperature of the isothermal object. Eq. (1) 
can be cast into a compact matrix form containing all nodes: 

rfV 

C — = AT + RT*" + A + R T;" + Q ( 2 ) 

dt S s s s 

where is a column vector containing the temperature at all nodes; is the temperature of the 

isothermal object. is the thermal capacitance matrix with the node -wise values C/ on its diagonal. 

A^^«x« is the conduction and radiation exchange matrices, respectively, filled with the 

values of the conductive and radiative conductors in the model. and are matrices 

storing the conductive and radiative links between the nodes and the isothermal objects (see Appendix for 
additional details). The superscript of *4 means the 4^^ power is applied to each element of a vector. 
Qg?{" is heat flux to each node. Since the heat flux Q can be split into constant flux Q^const and time 
varying flux Q^, we can define 

Qc=\^s+^sV+Qco.sr Q.=D« (3) 

i.e., Qc stores the constant thermal contribution from the isothermal objects and others (Qco^^/) on the 
nodes. Thus Eq. (2) can be rewritten as 

C— = AT + RT*^+Q^+Du (4) 

dt 

where u and D are, respectively, the independent inputs and the scatter matrix that converts the time- 
varying flux into thermal source to each node, such as the logic-controlled heaters and environmental heat 
flux. 


2.2 Reduced Order Modeling Algorithms 

Since the dimension of Eq. (2) is normally too large for efficient numerical integration, the ROM 
technique will be used to reduce the original full-scale model with a dimension ^ to a much lower 
dimension k by projecting the original model onto a low-dimensional subspace (i.e., T = UT^) 

constructed by a set of orthonormal basis vectors, where is the temperature in the reduced domain 

and k « n. Herein, we present two mathematically rigorous ROM algorithms, i.e., proper orthogonal 
decomposition (POD) and trajectory piecewise-linear (TPWL) reduced order modeling to identify an 
appropriate, low-dimensional projection subspace and automatically generate reduced thermal models 
amenable to fast computation by ODE/DAE solvers. 

Proper Orthogonal Decomposition Reduced order modeling (POD-ROM) 

Proper orthogonal decomposition (POD) [20] is a powerful method for data reduction and ROM, and has 
been extensively investigated for nonlinear dynamic systems, such as the rapid thermal processing 
systems [12, 13]. Therefore, it is an exceptional candidate to address strong nonlinear effects in spacecraft 
thermal analysis. For POD implementation, the method of snapshots is used, which extracts the subspace 
consisting of leading POD modes using an ensemble of data (snapshots) from the computation of the 
original, full-scale model [21]. To provide high-quality snapshot data, it is necessary to choose temporally 
characteristic, orbital input functions u (typically a series of sequences with relevant time scales) for the 
full-scale analysis. Given an input u, the matrix T containing the temperature snapshots at certain 
instances of time ti is given by 

r = [T(/,)T(»,)--T(/,)--T(v)] (5) 

where T(?/) is the /* snapshot in time; / = 1 ,2. . .A and N is the total number of the snapshots. Y then can be 
factorized using the singular value decomposition (SVD), followed by truncation at the first k « n 
leading modes yielding 
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where is the dominant subspace comprised of the k leading POD modes spanning the original 

snapshot matrix, = diag{c^} is a diagonal matrix of singular values arranged in a descending 

order. The singular value dj measures the dominance of its associated column in U. Therefore, by 
examining the relative magnitude of c^, an appropriate number k can be determined for ROM that yields 
desired model accuracy; and is the mode in the temporal domain. The original full-scale model 

Eq. (4) then can be projected onto the subspace U yielding a reduced thermal model with a much lower 
dimension (i.e., k « n): 

C,^ = A,T,+R^{UT,f+Q„,+D,u (7) 

where C, = U^CU, A, = U^AU, R, = U^R, = D, = U^D, T, = U^T . Given its 

significantly lower dimension, the ROM governing can be computed at much faster speed relative to 
the full-scale model using ODE/DAE solver. Figure 2 illustrates the flow chart of POD-ROM 
methodology. 
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C^ = f(T) + Du ^ C.-^ = U^f(UT ) + D u 


dt 


dr 


Figure 2. POD-based ROM algorithm development for spacecraft model. Snapshots from the full-scale 
model are factorized using singular value decomposition, followed by truncation at k leading modes, such 
that the ROM dimension k is significantly smaller than the full-scale model dimension n. 

Although POD is simple and accurate, its computational efficiency can be limited by the fact that the 
calculation of the nonlinear term requires converting the solution in the reduced domain back to the 
full domain T (as shown in Eq. (7)) and calculating the 4^^ order power of T on all nodes, leading to 
appreciable burden for computational acceleration (in particular for large models). To alleviate this issue, 
the discrete empirical interpolation method (DEIM) [22-24] can be used to approximate the nonlinear 
terms directly using the variable in the reduced domain, and hence, improves the computational 
efficiency. DEIM constructs specially selected interpolation indices that identify an interpolation-based 
projection to provide a nearly optimal subspace approximation to the nonlinear term. Therefore, only a 
few entries of the original nonlinear terms selected by the DEIM need to be evaluated at each time step. In 
DEIM, the nonlinear term is projected onto another subspace We with q « n, which in our case 
implies 

T*"«WH ^ (8) 

where W can be obtained using POD on the snapshot data of the nonlinear term only (i.e., 

(ti),...T*'‘(tN) used in Eq. (5) instead). H is the corresponding coefficient vector to be 
determined. Denote a matrix composed by q selected columns from the identity matrix 

i.e.,P = . If PW is nonsingular, the coefficient vector H can be determined from = 
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P^WH. Thus approximation to Eq. (8) becomes 

T*" » W (P^W)”' P^T*" = W (P^W)”' (p^t)*" . (9) 

The transformation in Eq. (9) allows the computation of the nonlinear terms only on a small number of 
selected nodes that have salient thermal contributions. Thus the reduced system of Eq. (7) can be 
modified as 

C, ^ = AX + R. W (p" w)“' (P"UT, )*' + Q„ + D^u = AX + X (PT, )*' + Q„ + D^u (1 0) 

where = R^W(P Wl , and P = P U can be pre-computed prior to the simulation. The 
interpolation indices Kp---,K^that construct the selection matrix P are selected inductively by the 
following algorithm [23]: 


Input; 

1. [/7Kj = max(|w,|) 

2. W = [w,],P = [e,J 

3. loop 1 = 2 to q 

Solve (P^W)z = P^w, forZ 

r = w,-WZ 
[p K,] = max(|r|) 

W = [Ww,],P = [Pe,J 

4. End loop 

Output: P for nonlinear function projection 

Figure 3. Algorithm for the discrete empirical interpolation method (DEIM) 

The function of “max” in the above algorithm returns the maximum component p in a vector along with 
the corresponding index of The process of the DEIM may be interpreted as selecting the indices of 
most representative components for approximating the nonlinear term [25]. 

Trajectory Piecewise-Linear Reduced order modeling (TPWL-ROM) 

Another effective means to address the nonlinear terms in ROM is to use the trajectory piecewise-linear 
(TPWL) technique [17, 18, 26-28]. That is, the full-scale nonlinear thermal model is first linearized along 
a series of linearization points. The TPWL approximation is utilized to construct a weighted combination 
of the linearized models to mimic the behavior of the original, full-scale nonlinear system in the entire 
domain. The linearized full-scale models in the weighted combination are then reduced to linearized 
ROMs using the subspace project method above. Specifically, the nonlinear thermal equation in Eq. (4) 
can be cast into a more compact form as 

irp 

C^ = f(T) + Du and f(X) = AT + RT*VQ^ (11) 

f(T) includes the conductive and radiative heat flux as well as the constant heat sources, and can be 
approximated using the Taylor expansion about a certain temperature T/, 

f(T)«f(T,) + h(T,)(T-T,) (12) 

where h(Ti) is the Jacobian of f(T) evaluated at T/ while T/ is the linearization point, i = 1, 2, ..., 5 . 
Note that substitution of Eq. (12) back into Eq. (11) yields a linearized full-scale model at the proximity 
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of T/. The original, full nonlinear model is then approximated by combining all the linearized models at 
the linearization points using TPWL method, i.e., 

C^=Z<><[f(T,) + li(T,)(T-T,) + Du] (13) 

at i=i 

where coi is the normalized temperature-dependent weight and = 1- Although linearized, direct 

computation of Eq. (13) is typically prohibitive due to its large dimension. Thereby a ROM procedure is 
normally taken to reduce its size, which is given as follows: determining a local projection subspace U/ at 
each linearization point in a trajectory/training simulation using a linear ROM technique (e.g., the Krylov- 
subspace ROM [29]). Then, an SVD is applied onto the union of all U/ to extract the global projection 
subspace where k again denotes the reduced dimension. Alternatively, the global projection 

space can also be obtained by directly applying POD/SVD on the snapshots obtained from the 
training/trajectory simulation. Based on our experience, the global subspace U obtained by the latter 
exhibits enhanced accuracy and robustness while the former allows subspace determination without 
resorting to full-scale model computation. Finally Eq. (13) is projected onto the global subspace U 
yielding: 


dt 


X <i(h, (T, ) T, + [f, (T, ) - h, (i; ) T, ] + ffl,D,,u 


(14) 


where fXT/)=U^f(T/), hXT/)=U^h(T/). Eq. (14) is essentially a weighted combination/interpolation of the 
linearized ROMs. The weights coi in Eq. (14) depend on the temperature of each node during the 
simulation, and can be calculated by the algorithm based on the Euclidean distance [17] as shown in 
Figure 4. Figure 5 illustrates the flow chart of the TPWL-ROM to generate reduced thermal models for 
spacecraft thermal analysis. 


1. For i= 1,...5 compute di = || J'-J'/IE. 

2. Take m = min {d^ for / = 1,. . .5. 

3 . For / = 1 , . . .5 compute W/ = 

4. Normalize w, set co^ = 

Figure 4. Algorithms to calculate the interpolation weights cOifor TPWL-ROM 
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Figure 5. Flow chart of TPWL-ROM 



Figure 6. Three different methods to locate the 
linearization points T/. 


An open question with TPWL-ROM is how to determine linearization points along a representative 
trajectory to accurately capture the dynamic response of the system. As illustrated in Figure 6, in general 
there are three methods to locate the linearization points. The first one {Constraint-\) uses a step-size 
constraint based on the simulation of the full-scale model, that is, the step size between two linearization 
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points (T/ and T/+i) needs to be below a tolerance l|T;+i- Till/ll Tjll < a^. In the second approach 
{Constraint-2), the step size is calculated using the linearized ROM and is restricted by 
||Tr,i+i “ ^ 2 - The third approach is essentially an error constraint {Constraint- 3), where 

the excursion from one linearization point to the next is determined in such a way that the error between 
the original full model and ROM is bounded, i.e., \\Tr,i “T^||/||TJ| < and hence, it requires the 
computation of both full-scale nonlinear and linearized ROM during the trajectory simulation. In this 
paper, the first and third constraints were utilized to locate the linearization points, and they typically took 
a value of 0.05-0.1. 

3. RESULTS AND DISCUSSION 

In this section, case studies using NASA-relevant spacecraft thermal analysis will be presented for ROM 
verification and demonstration. The original spacecraft thermal models in the form of the input data deck 
were exported from SINDA/FLUINT (S/F) [30], a widely used software tool in the spacecraft thermal 
analysis community. The S/F input data deck was then parsed using an in-house parser to extract pertinent 
model information and features, including: (1) nodal and inter-node information, such as capacitance and 
initial temperatures, conductive and radiative links/conductors, as well as the environmental and internal 
heating loads; (2) creating global indices for all nodes in various sub-models and for the nodal 
interactions; and (3) control logics used in the S/F model. The information parsed from the S/F model was 
then assembled in the form of Eq.(4), which, hence, is termed full-scale model hereafter. To quantitatively 
characterize the discrepancy between the full-scale model and ROM, two performance indices were used: 
the absolute error and the relative root mean square (rms) error, which are given by 

£"*=UT,-T and aT„,=7||UT,-T||/||T|| (15) 

where and T are the absolute temperature data from the ROM and the full-scale model. Note that at 
each time step Errabs is a vector and represents the temperature difference on all the nodes, while Err^ms is 
a scalar obtained by applying the L2-norm onto all the nodes in the model and indicates average relative 
error in the entire computational domain. To facilitate the comparison on the computational efficiency, 
both the full-scale model and ROM were simulated on the same platform, a Linux-based multi-user server 
equipped with 2 X AMD Opteron 6238 Twelve-Core 2.6GHz processor and 256 GB memory. The case 
studies using two spacecraft thermal models. Laser Interferometer Space Antenna (LISA) and Global 
Precipitation Measurement (GPM) with different computational sizes were analyzed. 

3.1 Laser Interferometer Space Antenna (LISA) Model 

The LISA mission, a space-based gravitational wave detector in the 0.1 to 1 mHz frequency band, 
measures distance fluctuations between proof masses aboard three spacecrafts. The LISA Integrated 
Modeling team developed a detailed thermal model that was used to drive the design of LISA [31, 32]. In 
this paper, a LISA thermal model containing 2874 nodes was used. In the transient analysis, the solar 
frequency is set at 0.1 mHz with 1% fluctuation in solar intensity (around the steady-state values). The 
overall heat source terms are given as 

Qc = Qeiec and Q, = Dmo (l + sin {27rcot)) 

where Qeiec is the electronics heating and Qt is the heat input from the solar radiation; uq and co are, 
respectively, the amplitude and frequency of the sinusoidal solar input. 

The model verification was undertaken in two steps: (1) comparing the S/F data from NASA against the 
full-scale model results to verify our in-house parser; and (2) comparing the full-scale model with ROM 
results to demonstrate salient computational efficiency achieved by ROM. Note that since the 
computational platform employed in (2) is different from that in NASA on which S/F analysis was 
performed, a direct comparison between S/F vs. ROM in terms of computational time and speed is not 
available. Figure 7 illustrates the transient temperature profile of all the nodes obtained using the full- 
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scale model and the S/F data from NASA with cd= O.l mHz to examine (1) above. The full-scale model 
exhibits excellent agreement with the S/F result in the entire spatiotemporal domain. The relative error 
Errrms in the entire temporal domain is <0.007% as shown in Figure 7c. The absolute error Errabs in the 
node -wise temperature and its temporal dependence is depicted in Figure 7d, which falls in the range of [- 
0.03K, 0.03K]. The temperature difference is negligible and may be caused by the different time stepping 
algorithms utilized by S/F and our full-scale model. Likewise, the transient simulation results with the 
solar frequency of co=\ mHz are very similar (not shown here). 
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(c) (d) 

Figure 7. Comparison between the full-scale model and S/F results in transient temperature profile for 
the LISA model with a solar fluctuation frequency of 0.1 mHz: (a) S/F data; (b) Data from the full-scale 
model; (c) Relative error; and (d) Absolute error in the node-wise temperature values. 


Next, the ROM is compared with the full-scale model in both accuracy and computational efficiency, i.e., 
step (2) above. In order to generate the snapshot data and ROM, trajectory simulation using appropriately 
selected training inputs that excite a broad spectrum of dynamic response in the nonlinear system were 
performed to determine the projection subspace U. The strategies to identify optimal inputs for snapshot 
simulation has been extensively discussed in prior research [33]. Normally it is an empirical process 
specific to individual systems, in particular, for complex environmental and heat source conditions (e.g., 
spacecraft). In our case study the training input is set as Wtrain ^ Uq/( 0? where Uo is the constant solar 
radiation intensity andXO is a time varying function [33] as shown in Figure 8. 
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Figure 8. Training input profile (ramp function) for solar radiation in the LISA model. 

POD-ROM 

The training input in Figure 8 was used in the full-scale simulation to generate numerical snapshots and 
extract projection subspace, which took 324 seconds in total. A projection subspace U with a low 
dimension of k= 12 and DEIM projection space W in rank ^ = 12 were then obtained for ROM. It should 
be noted that the training simulation and subspace generation is a one-time cost for reusable ROM, and 
hence, is negligible for the situation where ROM needs to be repetitively computed for parametric 
analysis. In addition, the training cost also heavily depends on the selection of the input profile. Although 
a simpler input with a short period can greatly shrink the computational time but at the cost of narrower 
operating envelop and limited reusability. In this case study ROM obtained using the input profile in 
Figure 8 is broadly applicable to all the simulation conditions of our interest (e.g., different solar 
frequency and varying amplitude). Figure 9a illustrates the transient, node -wise temperature profile, 
which matches very well the S/F data and the full-scale model shown in Figure 7a-b. Figure 9b-c, and 
Table 1 show the absolute Errabs and relative error Errrms of the ROM as compared with the full-scale 
model solution, which are less than [-0.6K 0.6K] and 0.02%, respectively. For computational efficiency, 
the transient analysis using the full-scale model cost about 266 seconds, and the POD-ROM achieved a 
speedup of 216x over the full-scale model. These findings confirm that the projection space obtained 
from POD and associated nonlinear ROM is able to accurately capture the strongly nonlinear heat transfer 
in the LISA model. 


Table 1. Computational performance for POD-ROM of the LISA model 


ROM dimension 

Errabs (K) 

Err yyyis /O 

Speed-up 

12 

[-0.6 0.6] 

<0.02 

216 



(a) (b) (c) 

Figure 9. Simulation results of POD-ROM (k= 12) vs. the full-scale LISA model: (a) Node-wise transient 
solution for ROM; (b) Relative error and (c) Absolute error. 
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TPWL-ROM 

Next we present the TPWL-ROM results and its computational performance for the LISA model. 
Different schemes to locate the linearization points along the trajectory in the training simulation were 
investigated. Figure 10 and Table 2 illustrate the comparison of the TPWL-ROM using constraint- 1 
(ROMi) and constraint-3 (ROM 3 ) against the full-scale model. The relative error Err^ms and absolute error 
Errabs of the ROMi are, respectively, <0.02% and in the range of [-0.4K 0.6K]. For ROM 3 , Err^ms is less 
than 0.06% and Errabs falls between [-0.4K 1.2K]. Note that the transient temperature profile obtained 
using the TPWL-ROM is almost the same as that by POD-ROM, and hence, is not shown herein for sake 
of brevity. In terms of computational efficiency, ROMi and ROM 3 exhibit a speedup of 186x and 213x, 
respectively, in contrast to the full-scale model. The ROM 3 runs slightly faster than ROMi due to fewer 
linearization points than the former although both have the same dimensions. 


Table 2. Computational performance for TPWL-ROM of the LISA model 



ROMi, a = 0.1 

ROM 3 , a = 0.05 

ROM dimension 

Errabs (K) 

Err yjyig /o 

Speed-up 

Errabs (K) 

Err yyyig /o 

Speed-up 

12 

[-0.5 2] 

<0.08 

186 

[-0.5 1.5] 

<0.12 

213 



Time(s) 

(a) 


x10" 



Time(s) 


Error AT(K) 



Time(s) 

(b) 


I 

I 


0.4 

0.3 

0.2 

0.1 

0 

0.1 

0.2 

0.3 

0.4 


xIO 


Error AT(K) 
1.2 



4 6 

Time(s) 


(c) (d) 

Figure 10. Simulation results of the TPWL-ROM (k=12) vs. the full-scale LISA model: (a, c) Relative 
error and absolute error, respectively, for ROM using Constraint-1 (ROMi); (b, d) Relative error and 
absolute error, respectively, for ROM using Constraint- 3 (ROM 3 ). 
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3.2 Global Precipitation Measurement (GPM) Model 

The next case study is the in-orbit analysis of the GPM satellite model obtained from NASA. The GPM 
model carries several distinct features and poses significant challenges to both full-scale model analysis 
and ROM. It has a large model size, containing 44,233 nodes and more than 1,700,000 conductors. Many 
nodes are made of composite materials and have temperature dependent properties. A total number of ~10 
orbits need to be simulated to eliminate the effect of the initial condition and examine its quasi steady- 
state thermal behavior. 

The GPM model has more than 100 heaters with logic control, and their input power depend on the 
programmed control law and time. Therefore, it is unrealistic to configure training inputs for each 
individual heat source to generate snapshot data. This is made even challenging due to the discontinuous 
and irregular profiles of the environmental flux as shown in Figure 11. Instead, we divided the orbital 
analysis into two stages: at the first stage, the first 1-2 orbits of the GPM model were simulated using the 
full-scale model and the snapshots data was collected and used to determine the projection subspace U 
and the low-dimension ROM as well. At the second stage, the rest of the orbits were analyzed using the 
ROM. The underlying rationale is that the temperature solution at the second stage can be treated as the 
perturbation to those at the first stage, and hence, the projection subspace will be almost the same for 
both. However, it should be pointed out that in contrast to the LISA model above, it will be difficult to 
reuse the ROM obtained by this approach for other orbit analysis that may have heat source profile 
notably different from the snapshot simulations. 

In the case study below, the first orbit of the full-scale GPM model was used as the snapshot data for 
subspace identification and ROM generation. The ROM was then used in the rest 9 orbits. In addition, our 
verification primarily focused on the comparison between extracted full-scale model and ROM with the 
S/F model neglected for two reasons: (i) the in-house parser has been verified extensively using the LISA 
and other relevant models; and (ii) the original GPM model in the S/F format from NASA included 
numerous heat pipes. However, HEATPIPE routines in S/F are not transparent to the external users, and 
cannot be parsed using our in-house tool at this stage. 



Timefs) 

Figure 11. An example of environmental flux applied to several nodes in the GPM model shows profiles 
with dramatic changes and sharp spikes. 

POD-ROM 

The ROMs with a dimension of ^ = 24 and ^ = 48 were first generated using the POD method. Figure 12 
portrays the transient node -wise temperature profiles obtained using the full-scale model and ROM with a 
dimension k = 24. The GPM model exhibits distinctly different thermal behavior at various regions as 
indicated by the dramatically oscillating temperature in some nodes while the rest remain almost constant 
due to strict thermal control requirements. As shown in Figure 13, the absolute Errabs and relative error 
Errrms is in the range of [-8K 8K] and less than 0.6% for the ROM with k = 24, respectively. However, as 
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the dimension increases to ^ = 48, the errors drop markedly, yielding -4K < ErVabs ^ 4K and Err^ms ^ 0.4%. 
It is also interesting to note that Errabs and relative error Errrms slightly grow with the time for both 
ROMs. It is caused by the shift in the subspace spanning the temperature solution between the initial and 
subsequent orbits. The simulation time of the full-scale GPM model in the ten orbits is about 80337 
seconds on our computational platform. The speedup of the ROM relative to the full-scale solutions is 
presented in two ways: the speedup was 9.8x and 9.2x if the computational cost of the snapshot 
simulation for the first orbit and the ROM generation was included. On the other hand, the speedup is 
104x - 458x in the second to the tenth orbit without taking this cost into account [34], which confirms 
that the ROM simulation time is negligible in contrast to the full-scale model. The quantitative 
performance metrics for the POD-ROM is summarized in Table 3. 


Table 3. Computational performance of POD in GPM model 


ROM 

dimension 

En-abs (K) 

En-ms (%) 

Speedup 

10 orbits 

0 
o 

1 

C3 

(N 

24 

[-8 8] 

<0.6 

9.8 

458 

48 

[-4 4] 

<0.4 

9.2 

104 



Figure 12. Node-wise transient solution of GPM model: (a) generated by full-scale model; (b) POD- 

ROM of dimension 24. 
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(b) 



Time(s) 


(C) (d) 

Figure 13. Simulation results of POD-ROM vs. the full-scale GPM model: (a, c) Relative error and 
absolute error, respectively, for the ROM with k = 24; and (b, d) Relative error and absolute error, 
respectively, for the ROM with k = 48. 


TPWL-ROM 

Table 4 and Figure 14 illustrate the comparison of the TPWL-ROM against the full-scale model solution 
using constraint-1 (ROMi) and constraint-3 (ROM 3 ). For k = 24, the relative error Err^ms and the absolute 
error Errabs are less than 0.8% and in the range of [-8 8 K]. As the dimension increases to k = 48, Errabs 
drops to [-5K 5K] and Errrms becomes less than 0.4%. In terms of computational efficiency, ROMi 
demonstrates a speedup of 9 . 6 x - 9 . 9 x with and 216x - 642x without considering the cost associated with 
the initial snapshot simulation and ROM generation cost as discussed above. Correspondingly, ROM 3 
exhibits a speedup of 9 . 5 x - 9 . 8 x with and 172x - 455 x without the cost associated with initial snapshot 
simulation and ROM generation. 


Table 4. Computational performance of TPWL-ROM for the GPM model 



ROMi, a = 0.1 

ROM 3 , a = 0.05 

ROM 

dimension 

(K) 

Erryyyis 

% 

Speed-up 

Errabs 

(K) 

Err^ffig 

% 

Speedup 

10 orbits 

2 ““ - 9* orbit 

10 orbits 

2 ““ - 9* orbit 

24 

[-8 8 ] 

< 0.8 

9.9 

642 

[-8 8] 

< 0.6 

9.8 

455 
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48 

[-5 5] 

<0.4 

9.6 

216 

[-5 5] 

<0.4 

9.5 

172 



Time(s) 


x10'^ 



Time(s) 



(C) (d) 

Figure 14. Simulation results of the TPWL-ROM (k=24) vs. the full-scale GPM model: (a, c) Relative 
error and absolute error, respectively, for ROM based on Constraint-1; and (b, d) Relative error and 
absolute error, respectively , for ROM based on Constraint-3. 


4. CONCLUSIONS 


In this paper, we presented a mathematically rigorous, subspace projection-based reduced order modeling 
(ROM) methodology for fast and efficient spacecraft thermal analysis. The detailed procedure of 
acquiring the whole spacecraft thermal model from the SINDA/FLUINT thermal analysis package and 
assembling it into an ODE/DAE form amenable for ROM was described. Two ROM methods, viz., 
proper orthogonal decomposition (POD) and trajectory piece-wise linear (TPWL), and their application to 
spacecraft thermal analysis for automated ROM generation were elucidated. In POD, the projection 
subspace is extracted from the solution snapshots using singular value decomposition (SVD) (or 
equivalent eigenvalue decomposition [35]), and the ROM is obtained by directly projecting the full-scale 
model onto the projection subspace. On the other hand, the TPWL-ROM is generated by weighted 
combination/interpolation of the localized linear ROMs to mimic the overall nonlinear dynamic behavior 
of the spacecraft. 

The computational performance of ROMs was demonstrated via case studies using NASA-relevant 
spacecraft thermal models, including the LISA and GPM. Depending on the nature of various thermal 
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models, different usage scenarios of the ROM were proposed. For the model with a small number of 
independent heat sources, e.g., LISA, carefully designed input function can be employed in the training 
simulation to reveal their individual contribution to the dynamic response of the system, and hence, the 
model generation in this case is a one-time cost and the generated ROM can be reused for various 
operating scenarios. For example, the LISA ROM can be used in various simulations involving different 
heat flux values (within the limit prescribed during the training simulation) and oscillating frequencies 
(0.1 and 1 mHz). The reusability of ROM enables several orders of magnitude acceleration in 
computational speed and significant resource usage relative to the full-scale model without appreciable 
compromise in simulation accuracy. For the model with a large number of independent heat sources, e.g., 
the GPM model, training simulation scanning through each individual source contribution becomes a 
formidable challenge. Therefore, spacecraft orbit analysis is targeted as the primary ROM application in 
this scenario, in which a small number of orbits (typically 1-2) are simulated using the full-scale model, 
the harvested data are used for subspace and ROM generation, and the ROM is then utilized for the rest 
orbits, typically yielding 5x - lOx speedup including the initial full-scale simulation ROM generation 
(>100x if these costs are not considered). In both case studies (LISA and GPM), POD-ROM and TPWL- 
ROM were verified by comparing against the high-fidelity, full-scale model, which exhibited excellent 
agreement in spatiotemporal thermal profiles (<1% relative error). 

The salient capability to capture spatiotemporal thermal behaviors that is otherwise unavailable by 
empirical models and the salient speedup over the full-blown analysis verified the utility of ROMs for 
initial concept screening, operational optimization, and thermal control and management. It should be 
noted that the present ROM is built on fixed operating and orbit parameters. To extend the current ROM 
framework to accommodate multi-parameter operations using sampling techniques for spacecraft design 
optimization [10, 11] will be the focus of our future work, which will effectively mitigate the 
aforementioned model reusability issues and significantly enhance the utility of the ROM. 
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Appendix: 

This appendix describes the process of assembling the node-wise information into A, R, A^, and 
matrices for thermal modeling of the whole spacecraft. and is filled with the values of 

the conductive and radiative conductors in the model as illustrated in Eq. (Al). Kij, 7?ij, are, respectively, 
the conductive and radiative links between the node and its Thermally’ adjacent nodes (via conduction 
or radiation). Ksi and Rsi are the conductive and radiative link between the node and an isothermal 
object (e.g., space), respectively. The value of Ksi and Rsi would be zero if the node has no flux exchange 
with the isothermal object. Therefore, A^ and R^ are diagonal matrices storing Ksi and Rsi for each node as 
Eq. (A2). 
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Subject to energy conservation, the entries in the row of matrix A and R need to satisfy 


(A2) 


n n 

T^4y=-Ki and 


(A3) 


This signifies that energy input to the node equals to that drawn from its ‘thermally’ adjacent nodes. 
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